/*
 * GeoDatabase.cpp
 *
 *  Created on: 13 Jul 2011
 *      Author: allan
 */

#include "GeoDatabase.h"

// C++ includes
#include <algorithm>
#include <iostream>
#include <fstream>
using namespace std;

// Boost includes
#include <boost/algorithm/string.hpp>
#include <boost/foreach.hpp>
using namespace boost;

// GeoReact includes
#include "Utils/StringUtils.h"
#include "Utils/Assert.h"

typedef vector<string> Block;

void ReadFile(vector<string>& strfile, ifstream& file);
void ReadBlocks(list<Block>& blocks, const vector<string>& strfile);

void ReadSpecies(SpeciesInfo& species, const Block& block);
void ReadSpeciesFormula(SpeciesInfo& species, const Block& block);
void ReadSpeciesName(SpeciesInfo& species, const Block& block);
void ReadSpeciesMolecularWeight(SpeciesInfo& species, const Block& block);
void ReadSpeciesVolumetricWeight(SpeciesInfo& species, const Block& block);
void ReadSpeciesRadius(SpeciesInfo& species, const Block& block);
void ReadSpeciesTags(SpeciesInfo& species, const Block& block);

void ReadComplexSpecies(list<ComplexInfo>& complexes);

void ReadReaction(ReactionInfo& rinfo, const Block& block);
void ReadReactionReactant(ReactionInfo& rinfo, const Block& block);
void ReadReactionProducts(ReactionInfo& rinfo, const Block& block);
void ReadReactionEquilibriumData(ReactionInfo& rinfo, const Block& block);

const string FindLine(const Block& block, const string& key);

const bool IsReactionBlock(const Block& block);

//-------------------------------------------------------------------------------------------------//

GeoDatabase::GeoDatabase()
{}

GeoDatabase::GeoDatabase(const string& filename)
{
	// Initialize the aqueous complex species
	ReadComplexSpecies(complexes);
	
	// Initialize the general species and the reactions
	ifstream file(filename);
	
	Assert(file.is_open(), "Database file could not be opened.");
	
	vector<string> strfile; ReadFile(strfile, file);
	
	list<Block> blocks; ReadBlocks(blocks, strfile);
	
	BOOST_FOREACH(const Block& block, blocks)
	{
		// Cheate a new species info
		SpeciesInfo sp; 
		
		// Read it from the current block
		ReadSpecies(sp, block);
		
		// Check if the current species in the block is an aqueous complex (if so, tag it with "complex")
		if(!FindComplexInfo(sp.formula))
			sp.tags.push_back("complex");
		
		// Add the new species
		species.push_back(sp);
		
		// Check if the current block is also a reaction block
		if(IsReactionBlock(block))
		{
			ReactionInfo reaction; ReadReaction(reaction, block); reactions.push_back(reaction);
		}
	}
}

const GeoDatabase::ReactionInfoList& GeoDatabase::GetReactionInfoList() const
{
	return reactions;
}

const GeoDatabase::SpeciesInfoList& GeoDatabase::GetSpeciesInfoList() const
{
	return species;
}

const GeoDatabase::ComplexInfoList& GeoDatabase::GetComplexInfoList() const
{
	return complexes;
}

const ReactionInfo GeoDatabase::GetReactionInfo(const string& reactant) const
{
	const optional<ReactionInfo> reactionInfo = FindReactionInfo(reactant);
	
	// Check if the information data of the species has been successfully found in the database
	if(!reactionInfo)
	{
		cerr << "Error: Cannot retrieve information of reaction whose reactant is " << reactant << ". There is no such reaction in the database." << endl; exit(0);
	}
	
	return *reactionInfo;
}

const SpeciesInfo GeoDatabase::GetSpeciesInfo(const string& formula) const
{
	const optional<SpeciesInfo> speciesInfo = FindSpeciesInfo(formula);
	
	// Check if the information data of the species has been successfully found in the database
	if(!speciesInfo)
	{
		cerr << "Error: Cannot retrieve information of species " << formula << ". There is no such species in the database." << endl; exit(0);
	}
	
	return *speciesInfo;
}

const ComplexInfo GeoDatabase::GetComplexInfo(const string& formula) const
{
	const optional<ComplexInfo> complexInfo = FindComplexInfo(formula);
	
	// Check if the information data of the complex has been successfully found in the database
	if(!complexInfo)
	{
		cerr << "Error: Cannot retrieve information of complex " << formula << ". There is no such complex in the database." << endl; exit(0);
	}
	
	return *complexInfo;
}

const optional<ReactionInfo> GeoDatabase::FindReactionInfo(const string& reactantFormula) const
{
	auto iterReactionInfo = std::find_if(reactions.begin(), reactions.end(), 
			[&](const ReactionInfo& reaction) { return reaction.reactant == reactantFormula; });
	
	if(iterReactionInfo != reactions.end())
		return optional<ReactionInfo>(*iterReactionInfo);
	else
		return optional<ReactionInfo>();
}

const optional<SpeciesInfo> GeoDatabase::FindSpeciesInfo(const string& formula) const
{
	auto iterSpeciesInfo = std::find_if(species.begin(), species.end(), 
			[&](const SpeciesInfo& species) { return species.formula == formula; });
	
	if(iterSpeciesInfo != species.end())
		return optional<SpeciesInfo>(*iterSpeciesInfo);
	else
		return optional<SpeciesInfo>();
}

const optional<ComplexInfo> GeoDatabase::FindComplexInfo(const string& formula) const
{
	auto iterComplexInfo = std::find_if(complexes.begin(), complexes.end(), 
			[&](const ComplexInfo& complex) { return complex.formula == formula; });
	
	if(iterComplexInfo != complexes.end())
		return optional<ComplexInfo>(*iterComplexInfo);
	else
		return optional<ComplexInfo>();
}

const GeoDatabase& GeoDatabase::GetInstance()
{
	static GeoDatabase db("Databases/chess.tdb");
	
	return db;
}

//-------------------------------------------------------------------------------------------------//

void ReadFile(vector<string>& strfile, ifstream& file)
{
	while(!file.eof())
	{
		string line; std::getline(file, line);
		
		//boost::erase_all(line, " ");
		
		boost::trim_if(line, boost::is_any_of(" \r\t"));
		
		if(boost::contains(line, "#") || line == "")
			continue;
		
		strfile.push_back(line);
	}
}

void ReadBlocks(list<Block>& blocks, const vector<string>& strfile)
{
	string current_species_type = "";
	
	unsigned deep_level = 0;
	
	BOOST_FOREACH(const string& line, strfile)
	{
		if(boost::contains(line, "{") && deep_level == 0)
		{
			current_species_type = boost::erase_all_copy(line, " {");
			
			++deep_level;
			
			continue;
		}
		
		if(boost::contains(line, "{") && deep_level == 1)
		{
			blocks.push_back(Block());
				
			blocks.back().push_back(current_species_type);
		
			++deep_level;
		}
		
		if(deep_level == 2)
			blocks.back().push_back(line);
		
		if(boost::contains(line, "}"))
			--deep_level;
	}
}

void ReadSpecies(SpeciesInfo& species, const Block& block)
{
	ReadSpeciesFormula(species, block);
	ReadSpeciesName(species, block);
	ReadSpeciesMolecularWeight(species, block);
	ReadSpeciesVolumetricWeight(species, block);
	ReadSpeciesRadius(species, block);
	ReadSpeciesTags(species, block);
}

void ReadSpeciesFormula(SpeciesInfo& species, const Block& block)
{
	species.formula = boost::erase_all_copy(block[1], " {");
}

void ReadSpeciesName(SpeciesInfo& species, const Block& block)
{
	string line = FindLine(block, "name");
	
	if(line == "") return;
	
	vector<string> words = Split(line, " =");
	
	species.name = words[1];
}

void ReadSpeciesMolecularWeight(SpeciesInfo& species, const Block& block)
{
	string line = FindLine(block, "molew.");
	
	if(line == "") return;
	
	vector<string> words = Split(line, " =");
	
	species.molecularWeight = atof(words[1].c_str());
}

void ReadSpeciesVolumetricWeight(SpeciesInfo& species, const Block& block)
{
	string line = FindLine(block, "vol.weight");
	
	if(line == "") return;
	
	vector<string> words = Split(line, " =");
	
	species.volumetricWeigth = atof(words[1].c_str());
}

void ReadSpeciesRadius(SpeciesInfo& species, const Block& block)
{
	string line = FindLine(block, "radius");
	
	if(line == "") return;
	
	vector<string> words = Split(line, " =");
	
	species.radius = atof(words[1].c_str());
}

void ReadSpeciesTags(SpeciesInfo& species, const Block& block)
{
	     if(block[0] == "basis-species")
		species.tags = {"aqueous", "basis"};
	else if(block[0] == "aqueous-species")
		species.tags = {"aqueous"};
	else if(block[0] == "gaseous-species")
		species.tags = {"gaseous"};
	else if(block[0] == "minerals")
		species.tags = {"mineral"};
	else if(block[0] == "redox-couples")
		species.tags = {"aqueous", "redox"};
	else if(block[0] == "surface-sites")
		species.tags = {"aqueous", "surface-site"};
	else
		cout << "Error: could not determine the type of the species " << species.name << "." << endl;
}

void ReadComplexSpecies(list<ComplexInfo>& complexes)
{
	ifstream file("Databases/complexes.dtb");
	
	Assert(file.is_open(), "Database file 'complexes.dtb' could not be opened.");
	
	string line;
	
	complexes.clear();
	
	while(!file.eof())
	{
		getline(file, line);
		
		boost::trim_if(line, boost::is_any_of("\n\t "));
		
		vector<string> words = Split(line, " ,=");
		
		ComplexInfo complex;
		
		complex.formula = words[0];
		
		for(unsigned i = 1; i < words.size(); i += 2)
		{
			complex.ions.push_back(words[i+1]);
			
			complex.stoichiometries.push_back(atof(words[i].c_str()));
		}
		
		complexes.push_back(complex);
	}
}

void ReadReaction(ReactionInfo& rinfo, const Block& block)
{
	ReadReactionReactant(rinfo, block);
	ReadReactionProducts(rinfo, block);
	ReadReactionEquilibriumData(rinfo, block);
}

void ReadReactionReactant(ReactionInfo& rinfo, const Block& block)
{
	vector<string> words = Split(block[1], " {");
	
	rinfo.reactant = words[0];
}

void ReadReactionProducts(ReactionInfo& rinfo, const Block& block)
{
	string line = FindLine(block, "composition");
	
	vector<string> words = Split(line, " =,");
	
	for(unsigned i = 1; i < words.size(); i += 2)
	{
		rinfo.stoichiometries.push_back( atof(words[i].c_str()) );
		rinfo.products.push_back(words[i + 1]);
	}
}

void ReadReactionEquilibriumData(ReactionInfo& rinfo, const Block& block)
{
	string line = FindLine(block, "logK");
	
	vector<string> words = Split(line, " =,()");
	
	for(unsigned i = 1; i < words.size(); i += 2)
	{
		// The log of the equilibrium constant log(K) (Note: the minus sign is necessary because CHESS write their reactions in the opposite direction we write
		double logKi = -atof(words[i].c_str());
		
		// The temperature at which the equilibrium constant is known
		double Ti = atof(words[i + 1].c_str());
		
		rinfo.logKi.push_back(logKi);
		rinfo.Ti.push_back(Ti);
	}
}

const string FindLine(const Block& block, const string& key)
{
	string line = "";
	
	for(unsigned i = 0; i < block.size(); ++i)
		if( boost::contains(block[i], key) )
			line = (boost::contains(block[i], "\\")) ? block[i] + block[i + 1] : block[i];
	
	boost::erase_all(line, "\\");
	
	return line;
}

const bool IsReactionBlock(const Block& block)
{
	string line = FindLine(block, "composition");
	
	return (line != "");
}

